! (c) 2008, Lawrence Livermore National Security Limited Liability Corporation.
! All rights reserved.
! 
! Redistribution and use in source and binary forms, with or without modification, are permitted 
! provided that the following conditions are met:
! 
!     * Redistributions of source code must retain the above copyright notice, this list 
!       of conditions and the following disclaimer.
!     * Redistributions in binary form must reproduce the above copyright notice, this list 
!       of conditions and the following disclaimer in the documentation and/or other materials 
!       provided with the distribution.
!     * Neither the name of the Lawrence Livermore National Security Limited Liability Corporation 
!       nor the names of its contributors may be used to endorse or promote products derived from 
!       this software without specific prior written permission.
! 
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR 
! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND 
! FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR 
! CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 
! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 
! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER 
! IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT 
! OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
      
      subroutine prec_scops(npoints,nlev,ncol,ls_p_rate,cv_p_rate, &
     &            cc,conv,cmip5_flag,precfrac_flag, &
     &            precfrac_ls,precfrac_cv,frac_out,prec_frac)


      use m_mrgzrnk

      implicit none

      ! Input variables
      INTEGER npoints       !  number of model points in the horizontal
      INTEGER nlev          !  number of model levels in column
      INTEGER ncol          !  number of subcolumns
      INTEGER precfrac_flag !  0 - no input large-scale or convective precip fraction
                            !  1 - use input large-scale precip fraction  
                            !  2 - use input convective precip fraction
                            !  3 - use input large-scale and convective precip fraction
      INTEGER cmip5_flag    !  0 - run the modified prec_scops
                            !  1 - run the prec_scops for CMIP5
      REAL ls_p_rate(npoints,nlev),cv_p_rate(npoints,nlev)  !  precipitation flux
      REAL cc(npoints,nlev) !  input cloud cover in each model level (fraction)
                            !  NOTE:  This is the HORIZONTAL area of each
                            !         grid box covered by clouds
      REAL conv(npoints,nlev)  !  input convective cloud cover in each model
                               !   level (fraction)
                               !  NOTE:  This is the HORIZONTAL area of each
                               !         grid box covered by convective clouds      
      REAL precfrac_ls(npoints,nlev)  !  input large-scale precip area fraction in each 
                                      !  model level
                                      !  NOTE:  This is the HORIZONTAL area of each
                                      !  grid box covered by large-scale precipitation
      REAL precfrac_cv(npoints,nlev)  !  input convective precip area fraction in each 
                                      !  model level
                                      !  NOTE:  This is the HORIZONTAL area of each
                                      !  grid box covered by convective precipitation

      REAL frac_out(npoints,ncol,nlev) ! boxes gridbox divided up into
                                       ! Equivalent of BOX in original version, but
                                       ! indexed by column then row, rather than
                                       ! by row then column
                                       !TOA to SURFACE!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      ! Output variable
      REAL prec_frac(npoints,ncol,nlev) ! 0 -> clear sky
                                        ! 1 -> LS precipitation
                                        ! 2 -> CONV precipitation
                                        ! 3 -> both
                                        !TOA to SURFACE!!!!!!!!!!!!!!!!!!
                    

      ! Local variables
      INTEGER i,j,ilev,ibox,cv_col,ib,thiscol
      INTEGER flag_ls, flag_cv  ! flag for completion of assignment
      INTEGER flag_p5_ls, flag_p5_cv  ! flag for possibility Five
      INTEGER sum_pf_cv
      INTEGER frac_out_ls(npoints,ncol),frac_out_cv(npoints,ncol) 
                       ! large-scale and convective cloudy levels in vertical columns
      INTEGER ind_ls(ncol),ind_cv(ncol)

      cv_col = 0.05*ncol ! assign default convective precip fraction (5%) into a number of subcolumns
      if (cv_col .eq. 0) cv_col=1

      if (cmip5_flag .eq. 1) precfrac_flag = 0
 
      ! initialize subcolumn precip fraction variable
      do ilev=1,nlev
        do ibox=1,ncol
          do j=1,npoints 
            prec_frac(j,ibox,ilev) = 0
          enddo
        enddo
      enddo
      
      ! count large-scale and convective cloudy levels in vertical columns
      do j=1,npoints
       do ibox=1,ncol
        frac_out_ls(j,ibox)=0
        frac_out_cv(j,ibox)=0
        do ilev=1,nlev
          if (frac_out(j,ibox,ilev) .eq. 1) then 
            frac_out_ls(j,ibox)=frac_out_ls(j,ibox)+1
          endif
          if (frac_out(j,ibox,ilev) .eq. 2) then 
            frac_out_cv(j,ibox)=frac_out_cv(j,ibox)+1
          endif
        enddo ! loop over nlev
       enddo  ! loop over ncol
      enddo ! loop over npoints

      ! working on the levels from top to surface    
       do j=1,npoints

        if ((precfrac_flag .eq. 1) .or. (precfrac_flag .eq. 3)) then
         call mrgzrnk(frac_out_ls(j,:),ind_ls)
        endif
        if ((precfrac_flag .eq. 2) .or. (precfrac_flag .eq. 3)) then
         call mrgzrnk(frac_out_cv(j,:),ind_cv)
        endif

        flag_p5_ls = 0
        flag_p5_cv = 0
        do ilev=1,nlev
        flag_ls=0
        flag_cv=0
     ! Large-scale precip assignment
        if (ls_p_rate(j,ilev) .gt. 0.) then 

          do ibox=1,ncol ! possibility ONE&TWO
            if (frac_out(j,ibox,ilev) .eq. 1) then !current level is cloudy
              prec_frac(j,ibox,ilev) = 1
              flag_ls=1
            endif
           if (ilev .ne. 1) then 
            if ((prec_frac(j,ibox,ilev-1) .eq. 1) &   !level above is precipitating
     &       .or. (prec_frac(j,ibox,ilev-1) .eq. 3)) then 
              prec_frac(j,ibox,ilev) = 1
              flag_ls=1
            endif
           endif
          enddo ! loop over ncol

          if ((flag_ls .eq. 0) .and. (ilev .lt. nlev)) then ! possibility THREE
            do ibox=1,ncol
              if (frac_out(j,ibox,ilev+1) .eq. 1) then  !level below is cloudy
                prec_frac(j,ibox,ilev) = 1
                flag_ls=1
              endif
            enddo ! loop over ncol
          endif

         if ((precfrac_flag .eq. 1) .or. (precfrac_flag .eq. 3)) then ! input ls precip fraction

          if ((flag_ls .eq. 0) .and. &
     &       (Int(precfrac_ls(j,ilev)*ncol+0.5) .gt. 0)) then
            thiscol=Int(precfrac_ls(j,ilev)*ncol+0.5)
            do ib=1,thiscol
              ibox=ind_ls(ib)
              prec_frac(j,ibox,ilev)=1
              flag_ls=1
            enddo
          endif
          ! if prec_scops precip fraction greater than model precip fraction, 
          ! clear columns, starting with those with fewest cloudy layers
          ! until fraction agrees
          if (Sum(prec_frac(j,:,ilev)) .gt. &
     &       Int(precfrac_ls(j,ilev)*ncol+0.5)) then 
          ! print *,'prec_scops overestimate'
            do ib=ncol,1,-1
              ibox=ind_ls(ib)
              if (Sum(prec_frac(j,:,ilev)) .eq. &
     &           Int(precfrac_ls(j,ilev)*ncol+0.5)) exit
              prec_frac(j,ibox,ilev)=0
            enddo
          endif
          ! if prec_scops precip fraction less than model precip fraction 
          ! (shouldn't happen often)
          ! start adding precip filled subcolumns until fraction is reached
          if (Sum(prec_frac(j,:,ilev)) .lt. &
     &       Int(precfrac_ls(j,ilev)*ncol+0.5)) then 
          ! print *, 'prec_scops underestimate'
            do ib=1,ncol
               ibox=ind_ls(ib)
               if (Sum(prec_frac(j,:,ilev)) .eq. &
     &            Int(precfrac_ls(j,ilev)*ncol+0.5)) exit
               prec_frac(j,ibox,ilev)=1
            enddo
          endif

         else ! no input ls precip fraction

          if (flag_ls .eq. 0) then ! possibility Four
            do ibox=1,ncol
              if (frac_out_ls(j,ibox) .gt. 0) then !large-scale clouds are in the vertical column
                prec_frac(j,ibox,ilev) = 1
                flag_ls=1
              endif
            enddo ! loop over ncol
          endif

          if (flag_ls .eq. 0) then ! possibility Five
           if (cmip5_flag .eq. 1) then 
             prec_frac(j,1:ncol,ilev) = 1  
           else
            if (flag_p5_ls .eq. 0) then  
             if((cc(j,ilev)-conv(j,ilev)) .eq. 0.) then
              if (((ilev .ne. nlev) .and. &
     &           ((cc(j,ilev+1)-conv(j,ilev+1)) .eq. 0.)) & ! if level below has no small cloud amount neglected by SCOPS
     &           .or. (ilev .eq. nlev)) then
               do ibox=1,ncol !large-scale precip covers 100%
                 prec_frac(j,ibox,ilev) = 1
               enddo ! loop over ncol
              endif
             else
               flag_p5_ls=1 ! precip flux with small cloud amount neglected by SCOPS
             endif
            endif  ! flag_p5_ls
           endif  ! cmip5_flag
          endif  ! flag_ls

         endif ! if large-scale precip fraction is input

        else
         flag_p5_ls=0
        endif  ! There is large scale precipitation
     
     ! Convective precip assignment
        if (cv_p_rate(j,ilev) .gt. 0.) then 

          sum_pf_cv=0

          do ibox=1,ncol ! possibility ONE&TWO
            if (frac_out(j,ibox,ilev) .eq. 2) then  !current level is cloudy
              if (prec_frac(j,ibox,ilev) .eq. 0) then
                prec_frac(j,ibox,ilev) = 2
              else
                prec_frac(j,ibox,ilev) = 3
              endif
              flag_cv=1
              sum_pf_cv = sum_pf_cv + 1
            endif
           if (ilev .ne. 1) then
            if ((prec_frac(j,ibox,ilev-1) .eq. 2)  &   !level above is precipitating
     &      .or. (prec_frac(j,ibox,ilev-1) .eq. 3)) then 
             if (prec_frac(j,ibox,ilev) .eq. 0) then
               prec_frac(j,ibox,ilev) = 2
               sum_pf_cv = sum_pf_cv + 1
             endif
             if (prec_frac(j,ibox,ilev) .eq. 1) then
               prec_frac(j,ibox,ilev) = 3
               sum_pf_cv = sum_pf_cv + 1
             endif
             flag_cv=1
            endif
           endif
          enddo ! loop over ncol

          if ((flag_cv .eq. 0) .and. (ilev .lt. nlev)) then ! possibility THREE
            do ibox=1,ncol
              if (frac_out(j,ibox,ilev+1) .eq. 2) then   !level below is cloudy
                if (prec_frac(j,ibox,ilev) .eq. 0) then
                  prec_frac(j,ibox,ilev) = 2
                else
                  prec_frac(j,ibox,ilev) = 3
                endif
                flag_cv=1
                sum_pf_cv = sum_pf_cv + 1
              endif
            enddo ! loop over ncol
          endif

         if ((precfrac_flag .eq. 2) .or. (precfrac_flag .eq. 3)) then ! input cv precip fraction

          if ((flag_cv .eq. 0) .and. &
     &       (Int(precfrac_cv(j,ilev)*ncol+0.5) .gt. 0)) then
            thiscol=Int(precfrac_cv(j,ilev)*ncol+0.5)
            do ib=1,thiscol
              ibox=ind_cv(ib)
              if (prec_frac(j,ibox,ilev) .eq. 0) then
                prec_frac(j,ibox,ilev)=2
              else
                prec_frac(j,ibox,ilev)=3
              endif
              flag_cv=1
              sum_pf_cv = sum_pf_cv + 1
            enddo
          endif
          ! if prec_scops precip fraction greater than model precip fraction, 
          ! clear columns, starting with those with fewest cloudy layers
          ! until fraction agrees
          if (sum_pf_cv .gt. &
     &       Int(precfrac_cv(j,ilev)*ncol+0.5)) then 
          ! print *,'prec_scops overestimate'
            do ib=ncol,1,-1
              ibox=ind_cv(ib)
              if (sum_pf_cv .eq. Int(precfrac_cv(j,ilev)*ncol+0.5)) exit
              if (prec_frac(j,ibox,ilev) .eq. 2) then
                 prec_frac(j,ibox,ilev)=0
                 sum_pf_cv = sum_pf_cv - 1
              endif
              if (prec_frac(j,ibox,ilev) .eq. 3) then
                 prec_frac(j,ibox,ilev)=1
                 sum_pf_cv = sum_pf_cv - 1
              endif
            enddo
          endif
          ! if prec_scops precip fraction less than model precip fraction 
          ! (shouldn't happen often)
          ! start adding precip filled subcolumns until fraction is reached
          if (sum_pf_cv .lt. &
     &       Int(precfrac_cv(j,ilev)*ncol+0.5)) then 
          ! print *, 'prec_scops underestimate'
            do ib=1,ncol
             ibox=ind_cv(ib)
             if (sum_pf_cv .eq. Int(precfrac_cv(j,ilev)*ncol+0.5)) exit
             if (prec_frac(j,ibox,ilev) .eq. 0) then
                prec_frac(j,ibox,ilev)=2
                sum_pf_cv = sum_pf_cv + 1
             endif
             if (prec_frac(j,ibox,ilev) .eq. 1) then
               prec_frac(j,ibox,ilev)=3
               sum_pf_cv = sum_pf_cv + 1
             endif
            enddo
          endif

         else ! no input cv precip fraction

          if (flag_cv .eq. 0) then ! possibility Four
            do ibox=1,ncol
              if (frac_out_cv(j,ibox) .gt. 0) then  !convective clouds are in the vertical column
                if (prec_frac(j,ibox,ilev) .eq. 0) then
                  prec_frac(j,ibox,ilev) = 2
                else
                  prec_frac(j,ibox,ilev) = 3
                endif
                flag_cv=1
              endif
            enddo ! loop over ncol
          endif

          if (flag_cv .eq. 0) then   ! possibility Five
           if (cmip5_flag .eq. 1) then 
            do ibox=1,cv_col
              if (prec_frac(j,ibox,ilev) .eq. 0) then
               prec_frac(j,ibox,ilev) = 2
              else
               prec_frac(j,ibox,ilev) = 3
              endif 
            enddo !loop over cv_col 
           else
            if (flag_p5_cv .eq. 0) then  
             if (conv(j,ilev) .eq. 0.) then
              if (((ilev .ne. nlev) .and. (conv(j,ilev+1) .eq. 0.)) &
     &           .or. (ilev .eq. nlev)) then
               do ibox=1,cv_col ! convective precip covers default fraction
                if (prec_frac(j,ibox,ilev) .eq. 0) then
                  prec_frac(j,ibox,ilev) = 2
                else
                  prec_frac(j,ibox,ilev) = 3
                endif 
               enddo ! loop over cv_col
              endif
             else
              flag_p5_cv=1
             endif
            endif  ! flag_p5_cv
           endif  ! cmip5_flag
          endif  ! flag_cv

         endif ! if convective precip fraction is input

        else
         flag_p5_cv=0
        endif  ! There is convective precipitation
        
        enddo ! loop over nlev
     ! end of initializing the top layer
       enddo ! loop over npoints

      end

